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Clothoids, i.e. curves Z(a) in R 1 whose curvatures kU) are linear fitting functions of arclength t, have 
been used for some time for curve fitting purposea in engineering applications. The first part of the paper 
deals with some basic interpolation problems for clothoids and studies the existence and uniqueness of 
their solutions. 

The second part discusses curve fitting problems for clothoidal splines, i.e. C 2 -curves, which are com- 
posed of finitely many clothoids. An iterative method is described for finding a clothoidal spline Z(») pass- 
ing through given points Z, t R 2 . i = 0,1 n+1, which minimizes the integral f K [>i 2 da. 

This algorithm is Buperlinearly convergent and needs only 0(n) operations per iteration. A similar 
algorithm is given for a related problem of smoothing by clothoidal splines. 

Key words: Approximation; clothoids; computer-aided design; Cornn-spirate; curvature; curve fitting; 
Fresnel-integrals; interpolation; splines 

Introduction 

The characteristic property of curves known as Cornu-spirals or clothoids is that their curvature x(s) 
is a linear function of the arc length, x(s) = k + As. Straight lines (x = 0, A. = 0) and circles (A — 0) 
may be considered as limiting cases. We are interested in constructing C 2 -curves in the plane R 2 which 
are composed of finitely many Comu-spirals; that is, C 2 -curves whose curvature is a continuous 
piecewise linear function of their arc lengths. We will call such curves clothoidal splines. Typical 
elementary problems encountered in such an effort are to construct a clothoid joining a given straight 
line and a given circle, or joining two circles. Composite curves of this type have been used by 
engineers, for instance, for the construction of highway sections, some of which are specified to be 
straight lines and circles. A more complex problem is to construct a clothoidal spline Z through a se- 
quence of finitely many points (at,-, v ; ) c R 2 , i = 0, 1 , . . . , n + 1 such that the integral 

K = fx[s) 2 ds 
Z 

along the curve is minimal. This problem can be considered as an approximation to the "true" prob- 
lem of curve fitting in it 2 , namely that of finding a curve Z(*) minimizing this integral among all 
C 2 -curves passing through the given points. The latter problem has been studied by several authors 
(Lee, Forsythe 17], 1 Mehlum [8]), and its exact solution leads to a multipoint boundary value problem 
for elliptic functions (Reinsch [14]). Mehlum [8] also proposed to approximate its solution by solving 
the corresponding multipoint boundary value problem for clothoidal spline functions, however the 
resulting clothoidal spline does in general not minimize the integral K among all interpolating 
clothoidal splines (see also Pal and Nutbourne [10] for a related use of clothoidal splines in computer 
aided geometric design). 

There is also the problem of smoothing: for given points (*,., y,), t = 0, 1, . . . , n + 1, the problem is 
to find a clothoidal spline Z in such a way that its deviation (in the least squares sense) from the given 
points is not greater than a prescribed tolerance and the integral K along Z is minimal (compare 
Reinsch [13] for the related problem for spline functions). 

*NBS Guest Worker with lie Operations Research Division, Center for Applied Mathematics, National Engineering laboratory. 
'Figures in brackets indicate literature references at the end of this paper. 
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Cornu-spirals can be easily computed in terms of Fresnel integrals, though admittedly not as easily 
as the cubic polynomials generally used for spline functions. In contrast to the latter, however, 
clothoidal splines are represented in terms of the natural parameter of plane curves; namely, the curva- 
ture as function of arc length. Furthermore, we hope that they do not exhibit the drawbacks observed 
with other schemes for curve fitting which have been observed in practice, namely, a tendency toward 
oscillations. 

In the first section we list some elementary properties of Cornu-spirals and Fresnel integrals, mainly 
taken from Abramowitz and Stegun [1], The second section deals with simple interpolation problems 
for a single Cornu-spriral. Section 3 is devoted to interpolation with clothoidal spirals; section 4 to the 
problem of smoothing. 

1 . Elementary properties of Cornu-spirals 

By definition, a Comu-spiral or clothoid is a curve, 



Z[s) = 



x{s) 
yis) 



,S£ R , 



whose curvature x\s) = x + As is a linear function of arc length s. If its tangent vector is 

. cos ${s) 



Z{s) = 



i <Ms 



then 



so that 



xis) = <M, 



Ms) = +0 + 1 tOOdr = +0 + *o» + ^ 



V 



(1.1) 



Z(s) = Z +f 



cos $(r) 
sin <{»(r) 



dt 



According to the sign of A, Z is called positively or negatively oriented. In the sequel, we restrict 
ourselves to the case of A > 0. Similar results will hold for k < 0. 
Using the Fresnel integrals, 



cos —— dt , S(z) := Jsin —r~ dt , F(z) := 
o 2 o 2 



CM := % 
o 

Z(s) can be expressed in closed form by [see [1], formulas (7.4.38), (7.4.39)] 

z W = Zo+ ^4-£)H^)- F (_^ 



F C\z? 
S(z) 



, if A > 0, 



(1.2) 



where V (a) is the orthogonal matrix, 
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tcos a -sin a~\ 
sin a cos a A 



Note that F[s) also describes a Gornu-spiral with arc length 5, curvature x(s) = its and phase angle <|>(s) 
= (n/2)s 2 . The Fresnel integrals have the following properties [see [1], (7.3.17), (7.3.20)] which we list 
without proof: 



Fiz) = -Fi-z) 



lim F(z) = - 

Z-»+oo 2 



, lim F(z) = - - 

Z—-co 2 



Moreover, F(z) can be expressed in the following way [see [1], (7.3.9), (7.3.10)]: 



m = \ 



-F(- Z 2)fc(z) , 



where the component functions g(z) and fiz) of 

hiz) = 



giz) 
fiz) 



satisfy [see [1], (7.3.5), (7.3.6), (7.3.21), (7.3.27M7.3.31)], 



(a) h(0)= - 



, lim h(z) = 



(b) g(z) and f{z) are strictly monotonically decreasing for z £ [0, + °°] 

(c) /' W = -msgW , g'(z) = nzf(z) - 1, for z e tf 

(d) For z > the following estimates hold for giz) and /(a): 



1 A. 




15 

(nz 2 ) 2 


-V giz) < ■ 


1 


*V ^ 


n 2 z 3 


1 ("■ 

ik y 




3 

(nz 2 ) 2 


-\< f(z)< 


1 

nz 


-3 

n 3 z° 


1 

JIZ 


.< 


TlV V 


35 


(nz 2 ) 2 



) 



(1.3) 



(1.4) 



(1.5) 



Approximations oif(z), giz) suitable for the calculation of F{z) are given in [1], (7.3.32), (7.3.33), and 
in Boersma [2]. 

As a simple consequence of (1.5d) we note the following estimates for the euclidean norms of the vec- 
tors h(z) and 



fib) := 
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giz) 

f{z) - l/inz) 



to be used later on: 



(a) 1- 



(b) 1- 



15 



(nz 2 ) 2 



< \\h(z)\\ (lAre) a/ 1 + 



f 



-1 



(us 



2(2 
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Xttz 



2)2 



< ||X(z)ll (1/irV Jl + 



nV y ; 



-1 



(to 2 )2 



<lforz>0, 



«lforz>0, 



<1.6) 



(c) lim zK(z) = 

Z-* + «> 



, lim m zh(z) = 

«— + « s --' 



By (1.5a), (1.5b), llftfzHI decreases strictly monotonically toward for z -»■ + °°, The same holds for 
H(z): 



I lK(z)l I decreases strictly monotonically toward as z -*■ + ». 
The following is a consequence of (1.6) and (1.5): 

lffifeJII 2 = —^m - — )<0 forz>0. 



2 dz 



ire" 



It follows from (1.3), (1.4) that F(z) has the form shown in figure 1. 



n 



0-5 




-0.5 ■ 



HKWll 




F(l) 



0.£ 



n »r 



FIGURE 1. Positively Oriented Comu-Spira! with Z = X = o and 



For z > the vector 

hiz)=\\Mz) 



cos o(z) 
sin o(z) 



>0, o(z) := arctan (flz)/«to) 
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(1.7) 



stays in the interior of the first quadrant otR 2 



0<o(z)< -.o(0) =~ , (+ oc> = n /2 

2 4 



Moreover, the vector 



r(z) := F(z) - 



"||- - v (p^ 



\\h(z)\\ 



cos o(z) 
sin q{z) 



where 



rotates counterclockwise for 2 > as z tends to + °°. This follows from (1.5): 

d fM 

q{z) = a[z) +nz = — arctant(/W/gfeW + *z = „. 2+^)2 



Therefore, the curve Fiz) crosses any fixed ray 

d 



>0forz>0 



■\l 





F 1 




r- - 1 




1 

2 


1 
1 


+ 


cob a 
sin a 





o>0 



infinitely often at abscissae < z\ < Z 2 < . . . , for which 

lim z; = + °° 
i-f 00 

(z i )2 + 4n-l«(z H . n > 2 <(z i )2 + 4n+l , i>l,n>l 

4n-l«(z n ) 2 «4n+l 

These estimates easily imply the following bounds 

4n+ 



tH 



1 + 



4n-l\- 1 



(1.8) 



< z i+n ~ x i < 



+1/ / 4^+iY 1 



V r 4n-l<z n < ^n+l , 

which we note for later reference. 

Upon inserting (1.4) into ( 1.2), we get the following representation of Z(s) in terms of the vector ft 



1 
(1.9) 



Z(s)= Z(0)- 



>[j(™*(£) -««*(£)) 



(1.10) 



where (see (1.1)) 



xW : = x + is, <Ks):= <t> + "0 s + ^ 



Note that because of A > and (1.5) (a), (1.3) 
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(a) Z(+~) =Z(0) + 



(b) Z(- oo) = Z(+ oo) 



-i\y(h- 



* 

2A 



(c) Z(*)-Z(+«) 



-^♦"•(w) 



(1.11) 



The evolute of Z, that is the locus of all centers of curvature M(s) of Z(s) for s £ R, is given by 



Mis)=Z[s) + 



x(s) 



-sin tys) 
cos $(s) 



= Z(») + 



xW 



nm 



(1.12) 



— VrM^)-^tS 



because 7i(z) = /i(z) ■ 



. Again, the evolute M is a spiral type of curve with the following properties: 



(a) M(+~) = Z(+~) 
<b> M(- »). = Z(- «) 



(c) MW-M(+«) 

(d) M(s 1 i*Mis 2 )ioTs 1 ^s 2 



= ~jl V{W {^ ] 






(1.13) 



(a) follows directly from (1.7) and (1.1 1), (b) and (c) follow from (1.11), and (d) from (1.7], since K(<|>(s)) 
is an orthogonal matrix. Furthermore, 



if x(s) > 0, A > 0, then for every s>s 
1 1 



IIAf(»)-AfWII < 



xi!) xis) 



(1.14) 



I IMG) - Z(s)ll < -— , 
x(s) 

that is, for s > s the osculating circle of Z at s and Z(s) are contained in the interior of the osculating cir- 
cle of Z at J. 

Indeed, according to a well-known result of differential geometry (see, e.g., [15]), the arclength ais) 
of the evolute Mis) of any curve Z(s) is given relative to the curvature x(s) of Z(s) by 



ois) = - —xis)- 1 
as 



so that in our case for s > s 
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o{s)-o(s) = — TZT 77- • 

xls) na) 

Since M(t), t e [s,s] is not a straight line, we have the additional inequality 

1 1Mb) -Af (s)l I < o(s) - o(s) = Kfi)- ! - icbr 1 
which proves the first part of ( 1 . 1 4) . The second part follows from the first, as 

NZW-Afla)ll = xW- x 
2. Interpolation properties of Cornu spirals 

In this section we study some simple interpolation problems for Cornu spirals. In stating the results 
we make use of oriented circles 



[r;]hH • 



Kia,r) :=\s + r . ,'||l<K+«2n 
whose orientation is determined by the sign of the radius r ¥= 0, and of oriented lines 

\oeR> , 



8 = 



ITcoso 
b + o . 
sin a 



whose orientation is deterined by the direction of the vector (cos a, sin a) ^*. We say that the orientations 
of an oriented hue g and of an oriented circle K(a,r) not meeting g are coherent, if K(a,r) lies in the 
same halfplane determined by g which contains the point 



b+r 




Coherent orientation 



Incoherent orientation 



Figure 2 
A first simple result refers to the problem of joining a line to a circle by a Cornu spiral. 

(2.1) THEOREM: 

1. For any given oriented circle K(a,r), r i= 0, not meeting a coherently oriented line g(b,a) there ex- 
ists exactly one oriented Cornu-spiral Z(s) which joins g to K(a,r) (in this order) such that the resulting 
composite curve is a C 2 curve with a coherent orientation, 

2. If g meets K or the orientation of g and K are not coherent, then there is no such interpolating 
Cornu-spiral. 

Of course, a similar result holds for joining an oriented circle K to an oriented line (in this order) by 
an oriented Cornu-spiral which we do not state explicitly. 
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PROOF: 1. Without loss of generality we may assume that r=l/x > and g is the x-axis in R 2 with its 
usual orientation. Since K(a,r) is coherently oriented with g, the center a = ix^y^) 7 of X is such that y 

;= y / r = y K > 1. 

Any positively oriented Cornu-spiral touching the x-axis at (0,0) r with s = (i.e., ^o = 0, Z(0) = 0) 
with a curvature x(0) = x = has the form (see (1.2)). 



Z(s) = 




V 




A 

-s I =: 



LrWJ 



with some A > 0. In order to solve the problem it suffices to determine s > and A. > such that Z has 
at s the curvature x and (x^y^ T as center of curvature (see fig. 3). 




This leads to the conditions 

x(s) = As = x 



Figure 3 



A = x/s , 



Us) =is 2 = >cs/2 , 
2 

cos 4>{s) = {y - y{s)) x = f-x x j '^- S 



m 



Hence s must satisfy the equation 



xs 



cos — + VnsJc S 
2 



W=)- • 



or the variable 



must solve 



Now the function 



Y : = 



cos 



v 2 + ^ V2T s (J- v ) = : 



p(^) : =.cobY 2 + 

= cos V 2 + 2VJ sin t 2 dt 


is strictly monotonically increasing for ^D because 

p' m = 2j sin Ait > for Y > 0. 


Since y~ > 1, p(0) = I and lim p(-r) = + °°, there exists therefore a unique solution V > of (2.2), 

T -mo 
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which can be found by Newton's method. In terms of V , the solution of the problem is 

s = 2¥ 2 /x , A = x/s 

72 



Z<s) =yjl F (4l3)>*o = *M-^j^ 



The proof of (2) is straightforward. 

We now turn to the problem of joining two oriented circles, 

K^a;, 1/k,), i= 1,2, , 
by an oriented Cornu spiral. 
We first show an auxiliary result for the family of Cornu spirals Z^ (s), A>0 with 

*o = <> : . 4o=°. z A (o) = o 



x(s) = As , $ is) = — s 2 



given by [see (1.10), (1.5a)] 



*»•- #("(?)»e)40) 



For their center of curvature M A (s) taken at arclength s: = x /A for which x(s) = x , the following 
holds: 



^-^HmMi) 



so that because of (1.6) (c), (1.11) and (1.13) 



(a) limM A (x/A) 

UO 



(b) limM A (x/A) + 

UO 



(c) lim M A (x/A) 
A-»+°° 



|_i/fJ. 



x>0 
if F<0 



(2.3) 



As an easy consequence, we get 



(2.4) THEOREM: Let K } ( ai , 1 /xi>, i - 1 ,2 fee two onented circles. 

1. J/Kj and K2 are coherently oriented, i.e. i/x; • X2 > 0, then there exists an oriented Cornu spiral 
joining Kj fo K2 (in this order) and having both K\ and K2 as osculating circle if and only if their 
centers aj are different and one of the circles contains the other in its interior. 

2, Ifxj • X2 <0, then there exists an oriented Cornu spiral joining Kj and K2 (in this order) and hav- 
ing both Kj and K2 as osculating circles if and only if neither circle contains the other, i.e. //aj -82!! > 
Kil|-7+ HK2II- 7 . 

PROOF: (1) Assume x^ > x\ > without loss of generality and let JCj contain K% in its interior; that is, 

0< II a!-a 2 II <1/jci-1/jc 2 ■ (2.5) 
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Then by (2.3) (a), (c) 

lim I! M A (X!/A)-M A (X2/A) II =0 
U0 

lim II M A ( X1 /A)-Mi(x 2 /A) II =l/x 1 -l/x 2 

A-~+« 

Since M k (x/A) depends continuously on A > 0, there is a A' > such that 

II M A ' (x!/A')-M A .()e 2 /A*) II = II aj-aj II , 

that is the Gornu spiral Z A , has two osculating circles of radii 1/xj and l/x 2 respectively, whose 
centers M A ' (xj/A' ), i = 1,2 have the desired distance. This proves the "if" part of (1). To prove the 
"only if" part, note that by (1.13)(d), the centers of curvature of any Cornu spiral are different for dif- 
ferent arclengths, so that aj + a 2 is a necessary condition for the existence of a Cornu spiral joining two 
different circles Ky, K 2 . The rest follows from (1.14). 

(2.) Assume Xj > > x 2 and II aj - a 2 II > 1/xj - l/x 2 . Then, because of (2.3) 

lim II M i (* 1 /A)-Af i (x 2 /J0 II = \/x l ~\/x 2 

A -. + 00 

lim II M A Ui/A)-Mi(x 2 /A) II = + » 

no 

Hence by a continuity argument there exists A' > such that 

II M k ' ix l /k')-M k 'ix 2 n') II ~ II a 1 -a 2 II 

which proves the "if" part of (2). The "only if" part is trivial. We next turn to the following problems: 

(2.6) PROBLEM: For a given oriented circle K and two points Pq £ K and Pj e K find an oriented 
Cornu-spiral connecting P to Pj (in this order) which has K as osculating circle at Pq (see figs. 4 (A), 
(B)). 

(2.6) is equivalent to the problem of connecting a point Py £ K to a point Pq t K (in this order} on an 
oriented circle K by an oriented Cornu spiral which has K as osculating circle at Pq. Using suitable 
reflections and changes of orientation [compare fig. 4 (B), (C)] ( (2.6) is seen to be equivalent to the 
following, which involves only positive orientations: 








** 



(A) (B) 




Figure 4 
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(2.6' ) PROBLEM: For a given positively oriented circle K = K(M , 1/x), k > 0, and two points Po £ 
K and P x £ K find a positively oriented Cornu-spiral with K as osculating circle at Po, which leads 
from P to Pi, if Pi is inside, K and leads from P\ to Po i/Pi. i/Pi is outside K. 

Clearly, (2.6' ) depends only on k and the relative positions of Pq and Pj so that we may assume 
without loss of generality 



M n = 



P* = 





-1/x 



>£l= ' 



sin a 
—cos a 



, r>0 , 



(see Fig. 5). 




Figure 5 



By (1.10) the class of positively oriented Cornu spirals Z with 

Z(0) 



is given by 



where 



= P =L 1/]d • ««» = * , +[o) = o 

= && (x/VS| - F(^( S ))A(x il (s)/"VS)) 

x A (s) : = x+As , ^(s) := xs + (A/2)s 2 . 



Essentially, we will show [Theorem (2.25)] that for r ^ 1/x, i.e. Pj c JC, there are countably many 
numbers Xj > X 2 > ■ •. • > and arclengths s,-, i > 1, such that Ci. («;) = P; for all i > 1. To prove this, 
we need some auxiliary results. From ( 1 . 5) (a) and (1.14) follow 



C k (+ ») 



-fr 



h ix/Tfnk) , II C,(+») II <l/x 



(2.8) 



We show next: 



(2.9) For any fixed bounded interval I = [sj, s 2 ] such that /or all A > and aK s £ I, x(s) = x + As >0 
there holds 



Hmsup II C A (s| 1 1 =1/k . 



(2.9) 
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PROOF. It follows from (2.7): 

C A ( 5 ) = J^ (h (x/ -VSX) - K(+ A W) A (-^)) " V^(s) 

By (1 .6)(c), the first two terms tend to uniformly in s e I as A 1 0. Hence, 



" 
l/x A W 



limsup IICjlWII = limsup l/jq(s) = 1/x ,QED. 



With the abbreviations 



h x :=J^= h (x/^d) , h x {s) :=J^- h(x x (s)/ nk) 
7 k := II h x II , r x \s) := II ft A W II , 



we have from (2.7) 



and from (1.6), (2.8), the estimates 

A /, 35A 2 



C k M=h-V[^[s))h x is) 



(2.10) 






< A/jc 3 , rl 



A 2 



F A < 1/x 

1 



for all s with jc a (») = x+ Xs > 0. 

Two cases are possible with respect to the location of the target point 

P: = r 



xAs) 



(2.11) 



sin a 
-cos a 



which will be treated somewhat differently. 

Case(l): 0< r< 1/x , Pj lies in the interior of K 
Case (2): r> 1/x , Pj lies outside of K. 

In Case (1) there is a sufficiently small A > such that 

C i (+oo)=h ;i ?tP 1 'forallO<A<A (2.12) 

Note this is exactly true if r = 0, P l = 0, for then by (2.8), 

C A (+ oo)=h i ^0 forallAX) . 
If r > 0, a suitable A > can be found because of (2.1 1). With T > satisfying (2.12), consider the rays 

d k -= {h k +oiPi-hJ\ o>0} , 0<A«A 
extending from h x towards Pj (see fig. 6). 
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W* }) 



Figure 6 



Because of (2.10) and using the same reasoning as with (1.8), every Cornu spiral C^s) , < A < A 
cuts d A infinitely often at abscissae < sj(A) <s 2 (A) < . . . , which satisfy estimates of the form [cf. 
(1.8)]. 



4> A (i n W) + 3n/2 < tikn+lUH 

2<n-l)n < 2(b- - fn < <Wa„W> < 2(n + - In < 2(n+l)n 

4 •*-*' 4 



for n > 1 , so that 



* fl+J (A)-* I1 (A) > 



3n 



x+As n (A) 



1 + 1 + 



J 



3nA 



[x+As„(A)] 2 



Vx<A) <*„(A» < ^ n+ i(A), 



where 



TJH : 



4m n 



1 + /1+. 



4mnA 



(2.14) 



(2.15) 



is the solution of the quadratic equation, 

4> A (s) = xs + — s 2 = 2mn . 
As a consequence, 

limsJA) = + oo, lim Cf (s„U)) = ftl > 

and therefore there exists an JV such that for all n > N (see fig. 6), 

Cx(s„(Al) £ [% , P x ) := {hX+ a (Pi-hj) I < a < 1 } , 

that is, Cjf intersects djj"between ft X and Pj at the abscissae s n (A ) , n < JV. Consider any fixed n < N, 
By (2.14), s„(A) is bounded 



m„ ^ s„(A) < M n for all < A < A 
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(2.16) 



by some positive constants m n , M n . Hence by (2.15), also the differences 

s„ +1 ii) - s„U) > m„ > f or all < A < I (2.17) 

are bounded below by a positive m„ > 0. Moreover, for each n>N, s n (A) is a continuous function of A, 
hence also C A (s„(A)), for 0< A < A . Since s n (A) is bounded above (2. 16), (2.9) gives for every fixed n 

wo 

that iB, the points P^ n := Q(s„(A)) £ d x tend to the boundary of the circle K as A tends to 0. Therefore, 
by the continuity of P x D and because of 

JPT.JS« lE.PJ 

there is a A„, < A„ < A such that P X t „ = Pj. 
Because of <2. 17), 

II C kn ( Sn+l (l n ))-h Xn II = a„K+ltt n >] < r in [s„(A„)] = II Pi-\H , 

so that 

Pi*Cx a ( Sn +i(W t ^ K ,Pi\ , 

and therefore A„ _|_ j < A„ . 

This proves that in case (1) there are indeed countably many positively oriented different Cornu- 
spirals C^ , n 3= 1, and abscissae s' n , namely 

having K as osculating circle at s = and passing through Pj, 

C x Js' n ) = P l , for all n 3=1. 

In caBe (2), r > 1/x, a similar reasoning applies: Here we consider the Cornu-spiral C x (s) for > s > 
-x/A, that is for all s < for which 

x k {s) = x + As > 

is still positive. We will show that: 

(2.18) To every integer n > 1 there exists a A > and an integer N > 1 such that for every < A < A 
the Cornu-spiral G^(s) , > s 3= - x/A, cuts d* at abscissae > s-i (A) > S-2U) . . . > 

(a) s_ N U) > -x/ F > -x/A , 

(b) s_j<A) - s^U) 3= m; > for i = 1,2 N+n-I, 0< A « X, (2.19) 

(c) r 7 [s_ N . 1 (X)]>r+- . 

(c) means that f or A = A , C x (s) has at least n cutting points, namely 
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Qrl*-jV-iU)l e {hI + a(Pi-hl)l a>l), i=l,2 n 

with dj" which lie beyond P\. 

Once (2J8) is proved, then as in case (1), a simple limiting argument MO gives the existence of n 
values A,-, A > A 2 > A 2 > . . . > A n > such that 

^(s-N-iik)) = Pi, 

since for A40 by (2.9) each C A . {j_jv— j(A», i > 1, tends to the circle Jf and so, by the continuity of s-jv-iU) 
has to pass the point P\ for a certain parameter value A;. 

Since by (2.18) n is arbitrary, this gives the existence of countably many Cornu-spirals satisfying the 
interpolation requirement. 

For the proof of (2.18) let y be defined by y/x := r + 1/x, so that y > 2. Let n > 1 be an arbitrary 
positive integer. Choose any numbers a and /J such that 



0<a<l , VT-fl«l/(2y) 
ajS < 1 , > 1 . 

Choose a natural number N so large that 

N+n+Kf}N 

a 2 * 1 



IVVll-al 2 2 

and set 



«2 



4Ntt 



Consider the solution s_ m , N j< m < fiN of the quadratic equation 



47" (j) = ks + — s 2 = -2mn 
A 2 



given by 



_ — imn {., , /. 4mnA i 
s -m= l* + J 1- — — ) 



x \ \j x 

Since by (2.20) 

every such s_ m is real. Moreover, 

A 7_jv = -ax/(\ + 1-a) = ~x ■ (1 - 1-a) 
A s _pjv = -x (1 - l-/}a-) 
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(2.20) 



(2.21) 



so that by (2.20) 

xjts-jv) = x + F-jv = k V r i^i > xtfs-pjv) = x \TFl)a > (2.23) 

Since by (2.21) 

15X 2 _ 15a 2 1 

xx(s-jv) 4 ~ 16iW(l-a) 2 2 

we get from (2.11) end (2.20) the estimate 

r X(s-jv) > r -i- t = -A >- = r +- . (2.24) 

x X(s-N) xyl—a x x 

Since by (2.21) 

k{s-p N ) - -2pNn < -2[N+ n+ 1 )ir , 

Cj{s) cuts d\ at least N+ n times within the interval [s-pjv, 0] at abscissae 

0>s- i a)>s- 2 a)>...>s- N - n (k) , 

satisfying the estimates 

-2ii-l)n>^(s-i(l)) >-2(i+l)n fori = 1,2, . . . ,N+n 

so that 

*_,*+ ] ^ s-,-(A) > s-i-i . 

In particular, we have > s_jv > 9 -jv-l(A), so that because of (2.24) and the monotonicity of rX(s), 
we get (2.19)(c). (2.19)(a) follows from s-N- n (X) > s-iv-n-1, (2.21) implying s-JV-n-1 > s-JJJV and 
(2.23). (2.19) (b) is proved as in case (1). All in all, we have shown the following: 

(2.25) THEOREM: For all oriented circles K and two points Pq £ K and Pi £ K there are countably 
many different Cornu-spirals connecting P to Pi (in this order) and all have K as osculating circle at 

3. Interpolation by Clotholdal Splines 

A clothoidal spline is a C 2 -curve in R 2 whose curvature x[s) is a continuous piecewise linear function 
of arclength s. More precisely, such a curve Z(s) is given by a finite collection of parameters 

= s <s 1 <. ..<s n+1 

(Z ; , <t> it xj, A,) , Zi£R 2 , 1 = 0, 1, , n 

such that for each i = 0,1, .... ,n , ZHs) : = Z(s)l[s,-, s,+ i] is a Cornu-apiral with curvature x'(s) and 
phase <r'(s) given by 
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ic'(s) 



:= X; + Xfis-Sj) 



*.-, 






(3.1) 



Z»'W : = 



s 

= z,. + ; 


cos 
sin 



(♦»'(*»* 



so that Z(s) is a c 2 -curve; that is, the Z l '(.), $'(.) , x'(.) satisfy the following continuity conditions for all 
i = 0, 1 , . . . , n — 1 : 



Z'( S( . +1 )-Z t .+ l s z f + J 



cos 
sin 



(<fr , '<S;+T})dT-Z < . + 1 =0 



A,- 



fVr !) - ♦,-+ ! s ♦,- + x f + x,- t,. + f*f- h+1 = ° 

X*k i+ 1) - X, + 1 = X; + A; T; - X,+ ] = 



(3.2) 



with t ( . := 8;+! — Sj. Of course, the parameters s,- are determined by the t,-, s,-_|_i = tq + tj + . . . + t,- 
so that instead of the *,-, we may take the t,- > as parameters. Note that we do not require A,- + 0, so 
that Z[$) may contain linear or circular segments. 

In this section we study the interpolation problem of finding a clothoidal spline passing through a 
finite number of given points. In this form, the problem is not very meaningful, since by Theorem 
(2.25) it has arbitrarily many different solutions. More interesting is the problem of finding an inter- 
polating clothoidal spline with minimal /x(s) 2 ds, in analogy to cubic spline interpolation. 

(3.3) PROBLEM: For a given family {Zj) i=0 (1 t n+1 of different points Z; E R 2 find parameters Pj 
= (<t>;, xj, Aj, tj), i = 0, 1, . . . , n with t ; > such that these parameters together with the Z; determine a 
clothoidal spline Z(s) by (3.1) satisfying (3.2) and Z(s n + 1 ) = Z n +i so that 

9 n+l n 3 i+l 

; x(s) 2 ds = z ; kWb 

i=0 s> 

is minimal. 
With the notation 



a?' : =^,x 1 ),6?':=(A,-,T 1 ) 
P T := (if, Pj, . . . , Pi) ,Pj:= (x t , Xj, A,-, t.) , i = 0,1 , 
the objective function to be minimized is the function 



(3.4) 



F(P):= I flxi+XppdT 
i=0 

which is separable in variables J*-. 

The transpose F' (P) of its gradient and its Hessian F"[P) are 
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F' (P) = (u , r , Ul , «!,..., u„, rj 



F"(i>) = 



(3.5) 



with the R 2 row vectors 



u,- := [0, 2x,t; + A,t?] 



V; := [Tf(x i +^A [ T I -)(K ( + A j T,)2] 



and the 4x4 square matrices 



(3.6) 



F r =2 



~0 
















T i 




-F 


, JC; + Xfii 





1 2 

' 2 T ' 




1 3 


, [Xf + AfTj-hf 





• "i + to 


ix 


+ V.K' 


, (X; + AfTjjA; 



(3.7) 



Also, the conditions (3.2) to be satisfied by P are highly structured. They have a staircase-like form 



G(P) = G{a 0t b a a ,b n ) = 



. K|a, b]) 



J(.„.,,h n . 1 H-Zi 1 -|-z B . 
K.(a n _„b n _il 



JlV^nl+Zn-Zn+l 
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where 



J(a,b) := / 




r 1 

cos 



sin 



j+ + *r + - T 2 WT( a:= 



V 


,b:= 


— - 


X 




T 



(3.9) 



K(a,b)-.= 



<|>+ XT+ - T 2 
it 

k + Xt 



Note that the Integral in J(a,b) is easily computed in terms of Fresnel integrals (see 1.2) for A; i= and 
elementary integration rules for A,- = 0. The Jacobian G ' of G has a similar structure 

G'(P).=G'(a ,b a n ,b n ) = 



^0 ' ^0 ' 
Co . Do , -J 

Ci.D!,-/ 







C n-1 , jDn-1 . -/ 



L 

with partial derivative 2x2 matrices 

Bi := JD^tt^THp. 



B n 







1 


» 


T f 


c i 


• — 















* 


1 



, Di := 



tf/2 , x { + AfT; 



In terms of the notation just introduced, (3.3) is equivalent to the minimization problem, 

Minimize F(P) subject to G(P) = 
Let 

L\P,A):=F\P)+A T G{P) 



(3.10) 



(3.11) 



(3.12) 
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be the Lagrangean of (3.12) and suppose that (3.12) satisfies the usual first order necessary arid second 
order sufficient conditions at the optimal point P (which we assume to exist): 

1. The Jacobian G' (P) of G at P has full row rank and there exists a A such that (P, A) is a sta- 
tionary point of L: 



A p L(P,A) 
G(P) 



L'IP.AF 



(3.13) 



<MP,A) = 0, with<r(P,A):= 

2. For the Hessian Lpp (P, A) of L with respect to P 

P T L pp (P,A)P>0 

holds for allF*0 satisfying G' (P)P = 0. 

Then P and A can be found as the solution of the nonlinear equations (3.13). The Jacobian +' of + is 
a highly structured matrix of the form 



+'(P,A) = 



L pp (P,A) , G'(P)T 
G MP) , 



(3.14) 



where G ' is given by (3.10). It is seen from (3.5), (3.10) that L pp has the same block-structure as F" 
(3.5). In solving (3.13), Newton's method can be applied to generate iterates (P' fc ', A' fc '), k = 0, 
1 , ... by solving at each iterate (P '*'• A " t ') the linear equations 



♦ ' (i* fc \ A<*>) 



r (JP<W' 
(JA*> 



tor the Newton direction 



tpm 

<5A<*> 



= -+(P*»,A*>) 



(3.15) 



, with <^' given by (3.14). 



Since computing the Hessian LpJP^ k \ A'*') may be too costly, we may replace Lpp within +' by a 
sufficiently close approximation JH™ as it is done in the minimization algorithms of Han [4, 5] and 
Powell [11, 12]. One may choose as H^\ e.g. a matrix of the same block structure as Lpp, namely 
(compare 3.5) 

~Hfr> 
Hf 



H»i=* 







H<*> 



(3.16) 



with 4 x 4 blocks iJ^' , i = 0, 1, . . . , n. One then solves (3. 15) with Lpp replaced by H^ k \ namely 

.n 



G'(P (ft ') 



G'iPW 




dA<W 



= H>(P ,fc) ,A<*>) 



(3.17) 
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and computes a new iterate of the form 

"p<fc+l)" 
\ik+l) 



p(fc) 
A [k) 



+ a k 



6 p\k) 



by choosing a step size a^ , 0< 04 < 1 , for example as in Han [5], by minimizing a certain penalty func- 
tion along the ray 



pik) 
A<W 



+ o 



6P m 
dA<w 



o>0 



After having computed the new iterate (p(*+D , A" t + 1 ') , one may use a rank-2 update formula, say 
the PSB-update formula, on each 4x4 block HW in order to generate another matrix H[ fc + " for each i 
= 0, 1, . . . , n, and thereby H* fc+1 ', having the same structure (3.16) as H**' and satisfying the usual 
Quasi-Newton equation: 



fltt+l)(pft+i» - Pf )) = V P Z(P< k + 1 U ,ft + 1> ) -V P .i(P (k, ,A (+1 >) 



(3.18) 



When solving (3.17), the structure of H' k ' (3.16) and G'(ffc) (3.10) can be exploited to reduce the 
number of operations drastically. For ease of notation, let us drop the superscripts and arguments in 
(3.17) and write briefly 

H 

for the right hand side -<f>(PM, A'*') of (3.17). The problem then is to solve an equation of the form 

(3.19) 



H , G' T 
G' , 



6P 
6\ 



r" 

c 



where H and G ' have the block structure (3.16) and (3.10), respectively. 

We first reduce G' by a series of Givens reflexions Qj, £ij* = Q,- , £ij = I, to a lower triangular 
matrix of the form [compare its structure with (3.10)]: 

G'- l .Q 2 ....Q N = (L,0) = 



\ 



\ 



O 




> 4nt2 



(3.20) 



where all blocks indicated have size 2x2 and L is a (4n + 2) x (4n -r- 2)-lower triangular band matrix. 
Again, because of the band-structure of (3.10), the number iV = O(re) of Givens reflexions needed is 
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linear in n, so that the unitary matrix 

Q := Qi'22 Qjv (3.21) 

need not be computed explicitly, but can be stored in product form. Partition the matrix 

Q = (Q, £2) 
where 



Q=Q 



are the last two columns of Q, which are computed using the product form of (3.21), Q is not needed 
explicitly. Introduce new variables 













1 
1 


= Qj.Qa... 


..Qjv 


1 
1 



t = 



via 6P = Qt = Qfj + Qt 2 . 
Then because of 

G'Q = L,G'Q=0 
the second set of equations (3.19) 

G'6P = Lt l = d-"t 1 
can be solved for t j in 0(n) steps using the structure of L (3.20), and the vector 



P 1 := Qt Y = Q X Q 2 . . . Q N 

is computed using (3.21). 
Now we turn to the first set of equations (3.19) 

H6P+G' T 6\ = c 



h 





(3.22) 



Multiplying these equations by Q T and introducing tj and t 2 instead of 6P, we get because of 

QFHQti + QTHQt 2 = Q T c 



or 



{Q T HQ)t 2 = Q T c - Q T HPl -~ 1 2 
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(3.23) 



Again, the 2x2 matrix Q T HQ and the vectors Q T HP^ can be computed with 0(n) operations using the 
block structure of H (3.16). r 2 is obtained by solving the two linear equations (3.23) and 6P is 
calculated by 

P 2 :=Qt 2 ,6P:=P 1 + P 2 . 

Finally, we multiply (3.22) by Q T in order to get <JA. Observing (3.20) we obtain a triangular system of 
linear equations 

L r t<JA = Q r c-2 r ffdP 

the right hand side of which can be easily computed with 0(n) operations using the structure of H and 
the product form of £2^ 



2T = 



1 







000 



100 



QNQjv-j . . . Q 1 



All in all, we can compute the solution of (3.19) with 0(n) arithmetic operations, so that the Han- 
Powell method is quite effective in our case. The method has been realized and successively tested by 
Huckle [6]. With respect to a convergence analysis of the above method (the method converges locally 
superlinearly under some mild assumptions) we refer to the literature Han [4,5], Powell [12], Tapia 
[16]. 

4. Smoothing by Clothoidal Splines 

We consider the following generalization of (3.3) (compare Reinsch [13]): 
PROBLEM: For a given family ^L^ i=LfiX n+1 of different points 

£R2 



Zi = 



7i 



and numbers S > 0, Axj > 0, Ayi > 0, i = 0,1 , . . . , n+ 1 , find parameters 



{{(<l>i,Xi,Ai,Ti))i=0,l D »{ Z i}i=0,l n + 1, z}, Z; = 



:R2 



which determine a clothoidal spline Z(s) via (3.1) satisfying the conditions 
a) (3.2)andZ(s n+1 ) = Z n+1 



(4.1) 



(4.2) 



" tm + m 



+ z 2 = S 



[zisa 



slack variable) such that J A(s) 2 ds is minimal. 
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Again with the notation [compare (3.4)] 

P„+i T = [Zo«,Z 1 r... f Z Il+1 T, z ] 

= (x .yo»*l.yi x n+1 ,y n+1 ,z)ER 2n+5 

the objective function F(f1 to be minimized is separable in the P; 



i= 0,1, . . . , n 



F(P): = Z / (x, + A,t) 2 dr 

i=0 



and has a Hessian of the form [compare (3.5)] 



F m iP) = 



•F„ 







'n+l 



(4.3) 



(4.4) 



with the same 4 x 4 matrices F F„asin(3.7)anda(2n+5)by (2n+ 5) matrix F n+ i 

The constraints (4.2) now have the structure [see (3.8)]: 

J(aQ,b$) + Z - Z\ 

J(ai,6l)+Zi-Z 2 
M*l,bi)-a2 



= 0. 



G|P) = 



JK-ubn-rf+Zn-l-Zn 
K K-l,b n -i)-a n 
JK,b n )+Z„-Z n+1 
C(Z ,Zi„..^Z n+1 ^) 

where L and K are again given by (3.9) and the scalar function q is defined by [compare (4.2) b)] 



= 



(4.5) 
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« z ° ^ =p^f#- 

With these new definitions of F and G, problem (4. 1) has the same structure as (3. 1 2), namely 

minimize F[P) subject to G(P) = (4.6) 

Consider again the Lagrangean of (4.6) 

UP,A):=F[P) + A T G(P) 

We again assume that (4.6) has an optimal solution P and that at P the optimality conditions (3.13) are 
satisfied. 



By (3.13), (3.14) the optimal solution (P, A) solves 

*(P,A):=L'(P,A) = 
whose Jacobian is again 

f (P,A> = 



ApL(P,A) 
G{P) 



= 



(4.7) 



Lpp(P,A) , G'T 

G' , 



(4.8) 



but its structure is slightly more complicated than in section 3 because of our new definitions of F[P) 
(4.3) and G(P) (4.5). It is easily verified that in the present case $'(P,A) has the following form 
(illustrated for n=2) 



} 2 



G'(P) 



f— 








" 


Ao,B 


I, 


-I, o. 








Cn.Do. -1 


0, 


0, 0, 








Aj.Bj 


0, 


I, -I, 








Ci.Dj, -I 


0, 


o, o, 








A 2 ,B 2 


0, 


0, I , 


-I 








r 


z 



(4.9) 



} 2 
} 1 



2 2 



where the 2 by 2 matrices A it Bi, Q, D; are again given by (3. 1 1 ) and the vector r is 



r:= 



*o~*b yo - 7o 



A* > 2 (Ayo» 2 



»n + l-g/i+l 
(A*„+l> 2 



yn+i~y n +i 

(Ay n+ ,) 2 



Likewise L PP (P,A) has the structure [compare (4.3), (4.4)] 
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L PP (P,A) = 



(4.10) 







J n+1 



where the L it i < a, are symmetric 4 by 4 matrices and L n +] is the (2/i+5) by (2n+5) diagonal 
matrix. 

where A z is the last component of A. 

As in the previous section, one has to solve (4.8) by Newton's method (compare (3.15) - (3.17) 
where at each iteration point [P" 1 *, A" 1 '] the Hessian L PP is approximated by a positive definite matrix 
H' fe ' having the same structure as Lpp (4. 10), 



H ik): 



H W 



jj.tft) 



rV"» 







H n+ i< fc > 



(4.12) 



with certain 4 by 4 matrices H/&) f or i < n and the diagonal matrix (see 4. 1 1 ) 
H n+J M = Aj.< fe > * diag(Ax ,Ayo, .... Ax„+i,Ay„+i,l)" 2 



(4.13) 



After having computed P^ h+l \ A ,fe+1) (see previous section) H <fe+1 ' is obtained from H ,fe> by updating 
each H* k \ i < n, individually by some update method (e.g., the PSB-method) which guarantees the 
same quasi-Newton relation (3.18) as m section 4; H^^.^ k+l ^ is computed by (4.13). 

Of course, for large numbers n the efficiency of the algorithm outlined crucially depends on the 
number of operations needed to perform one Newton step [P ,fc) , A< fc) ] -* [P Wc+I> , A (fe+1 >], that is to solve 
a linear system of equations [see (3.17), (3.19)] of the form 



G' ,0 



6P 

dA 



(4.14) 



for dP, dA, where H and G' are given matrices with the structure (4.12) and (4.9), respectively. An 
algorithm of the type considered at the end of the previous section leads to difficulties inasmuch as it 
would take 0(n 3 ) operations to solve (4.14) because it requires the computation and storage of a large 
dense matrix of the order 0(n). 

Another numerically stable way to solve the linear system (4. 14), which exploits the symmetry of the 
matrix 



H , G' T 

G' , 
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(4.15) 



would be to use the Bunch-Parlett decomposition of (4.15) (see Bunch, Parlett [4]). However, this 
method requires a pivot selection in each basic elimination step, which, though preserving the sym- 
metry, will in general destroy the specific block structure of the matrix in (4.15). This method, 
therefore, also requires 0(n 3 ) operations to solve (4.14). A cheaper method for solving (4.14) might be a 
variant of the conjugate gradient algorithm for solving linear equations 

Ax = b 

with a symmetric nonsingular, but perhaps indefinite matrix A, which is described in Paige and 
Saunders [9]. This method can take the block structure (4.12), (4.9) of H and G' into account and 
therefore requires only 0(n 2 ) operations and 0(n) storage to solve (4. 14). 

It is interesting to note in this context that the system (4.14) can be solved with only 0(n) operations, 
if the block-diagonal matrix L PP (p, A) (4.10) would be positive definite at the solution (P, A) of (4.7). In 
this case, it can be shown that the matrices iJ* fe> (4. 12) generated by the usual update techniques (PSP-, 
DFP-, or BFGS-methods) wiU be positive definite, at least locally, if the starting values [P* *, A (0> , and 
H* 0> ] are sufficiently close to (P, A) and L PP (p, A), respectively. 

If H is positive definite, then a numerically stable method of solving (4.14) requiring only 0(n) 
operations runs as follows: 

In a first step compute the Cholesky decomposition of 

H = RTR 

which requires 0(n) operations and gives an upper triangular R of the form [compare (4.12)] 



R = 



-Rn 







it, 



R. 



(4.16) 



with 4x4 upper triangular it,- for i < n and diagonal R n+ 1 . Premultiplying (4. 1 4) by 

> T , 0~ 

-G'R- l R- T , I 



gives the equivalent system 

~R , (G-R- l ) T 



6P 
6/\ 



ii-Tc 
d-G'R^R-Tc 



(4.17) 



-{G-R-mG-R- 1 ) 7 

So the next step is to compute 

c'-^R-TcA^G-R-l 

which again requires only 0(n) operations because of the simple structure of R (4.16) and G' (4.9). 
Note, moreover, that the product matrix A — G 'it -1 has a form very similar to (4.9), namely (illustra- 
tion for n— 2): 



343 



x x x x 

X X X X 

X X X X X o 

X X X X X X 

X X X X 
X X X X 
X X X X X o 
X X X X X X 

X X X X 

x x x x 


x o x o 

O X O X 

X X o 

O X O X 

X X 

o x o x 


-I 



o 
o 





xxxxxxxx 


x 



(4.18) 



We next reduce A to "lower triangular" form by multiplying A from the right by suitable Givens 
reflexions Qj, Q 2 > • • • > &N = ^M matrices 2,- and only 0(n) operations are needed and the structure of 
(4.16) is essentially preserved and fill in will occur at most 0(n) places. Each will annihilate a particular 
above diagonal element of A; the resulting matrix is of the form 



AQ&2 . . . Qn = (L, 0) 



(4.19) 



where "0" denotes a (4n+3) by (2n+6) zero matrix and L is a (4n+3) by (4n+3) lower triangular 
matrix with the structure 



l - 




[>' 



■. k 



ft 



^ 



Note that the dense (6n+9) by (6n+9) product matrix 2=2x22 . • • • » 2 Af Q eed not be computed. 
Its storage in product form requires only 0(n) places. Concurrently with the elimination process for 
finding L, we can compute the vector 

c":=2 Jv ...2 2 2 1 c' 

Now it is easy to solve (4. 1 7) for iP and dA. The second equation (4. 1 7) gives by (4. 1 9) at once 

AA T 6\ = LL T 6A = -d + Ac' = -d + <L,0)c" (4.20) 

so that 

dA = -L-TL-id + {L- T ,0)c" (4.21) 

i.e. dA can be found by solving three linear equations with triangular matrices. The first equation(4. 1 7) 
now gives by (4. 2 1 ) 
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R6P= R~ T c~A T 6A 




= c 



tJA 



(4.22) 



= c'-Q 

Unfortunately enough, the computation of 

L-ld-(/,0)c" 



:=Q 







~L- l d+iI,0)c' 




= QjQj , . . . , Qjv 



r L- 1 d-(/,0)c* 




requires the storage of all £2,- (this was not needed in computing c"). Note that L~^d has already been 
obtained during the calculation of (JA (4.21). Finally, by (4.22), 6P is obtained by solving one more 
triangular system of linear equations 



R6P = c , -c'"^6P , 



(4.23) 



again requiring only 0(n) operations. 



At the expense of numerical stability one may get around the elimination process to find L and the 
storage of the orthogonal matrices Q,- in the following way: 

Having computed the Cholesky decomposition of H = R T R, the matrix A = G 'it -1 , the product 
AA T and its Cholesky decomposition AA T = LL T , computing <JA and 6P from (4.20), (4.22) is 
straightforward: 



LL T 6A = -d+Ac'-6A->A T 6A 
R6P = c'~A T 6A^6P 



(4.24) 



Note in this context that the product AA T has a simple sparse structure needing only 0(n) places for 
storage: 



AA* 
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Both algorithms require only 0(n) operations for solving (4.14) in each Newton step, but the former 
will be numerically more stable, as it avoids the calculation of AA^- and cancels products such as 
LLT 1 , RR -i , which arise inherently during the solution of (4.24), as often as possible. 
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